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Abstract 

A  non-isothermal  dynamic  optimization  model  of  direct  methanol  fuel  cells  (DMFCs)  is  developed  to  predict  their  performance  with  an  effective 
optimum-operating  strategy.  After  investigating  the  sensitivities  of  the  transient  behaviour  (the  outlet  temperature,  crossovers  of  methanol  and 
water,  and  cell  voltage)  to  operating  conditions  (the  inlet  flow  rates  into  anode  and  cathode  compartments,  and  feed  concentration)  through  dynamic 
simulations,  we  find  that  anode  feed  concentration  has  a  significantly  larger  impact  on  methanol  crossover,  temperature,  and  cell  voltage  than  the 
anode  and  cathode  flow  rates.  Also,  optimum  transient  conditions  to  satisfy  the  desired  fuel  efficiency  are  obtained  by  dynamic  optimization.  In 
the  developed  model,  the  significant  influence  of  temperature  on  DMFC  behaviour  is  described  in  detail  with  successful  estimation  of  its  model 
parameters. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Direct  methanol  fuel  cell  (DMFC)  systems  are  attractive  for 
portable  and  stationary  power  applications.  The  fuel  (methanol) 
is  easy  to  store  and  deliver,  and  has  a  high  specific  energy.  The 
design  complexity  of  DMFC  systems  is  reduced  by  eliminating 
the  reforming  units  that  are  required  to  produce  hydrogen  for 
polymer  electrolyte  membrane  fuel  cells  (PEMFCs).  In  view  of 
designing  the  systems,  there  are  two  types  of  DMFCs,  namely, 
passive-feed  DMFCs  and  active  DMFCs.  Studies  of  passive 
DMFCs  have  been  undertaken  by  Chen  and  Yang  [1],  Kim  et 
al.  [2],  Liu  et  al.  [3],  Shimizu  et  al.  [4],  and  Chen  and  Zhao 
[5].  Passive-feed  DMFCs  have  advantages  in  eliminating  para¬ 
sitic  power  loss  of  pumps  and  fans  because  without  an  active 
(convective)  feed  supply,  oxygen  is  supplied  by  air-breathing 
and  methanol  is  diffused  into  the  anode  from  the  fuel  cartridge 
driven  by  a  concentration  gradient.  By  contrast,  active  DMFCs 
can  control  the  flow  rate  of  the  feed  and  air  by  means  of  pumps 
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and/or  fans  to  cope  with  the  varying  power  loads  and  environ¬ 
mental  conditions.  However,  both  types  of  DMFCs  still  have  low 
fuel  efficiency  and  specific  power  compared  with  hydrogen  fed 
PEMFCs  due  to  mass  crossover  problems  and  slow  methanol 
oxidation. 

The  fuel  efficiency  and  the  available  oxygen  concentration 
at  the  cathode  catalyst  sites  are  reduced  by  methanol  crossover 
through  which  a  large  amount  of  methanol  is  transported  from 
the  anode  to  the  cathode  catalyst  layer.  Oxygen  transport  to  the 
catalyst  layer  is  also  prevented  as  water  at  the  anode  crosses 
the  membrane.  To  solve  these  problems,  much  research  has 
been  conducted.  In  modelling  work,  Simoglou  et  al.  [6,7]  have 
described  an  experimental  study  of  the  dynamics  of  a  DMFC  and 
have  developed  an  empirical  model  of  the  cell  voltage  dynamics 
with  a  feasibility  study  of  model-based  scale-up  and  scale-down. 

Dohle  et  al.  [8]  examined  the  heat  and  power  management  of 
a  DMFC  system,  and  in  particular  investigated  the  influence  of 
water  vaporization  in  the  cathode  on  the  heat  management  of  the 
DMFC  system.  Dohle  and  Wippermann  [9]  developed  a  semi- 
empirical  model  to  predict  methanol  permeation,  overpotentials, 
and  polarization  curves.  Sundmacher  et  al.  [10]  formulated  a 
dynamic  lumped  parameter  DMFC  model  by  deriving  anode 
kinetics  consisting  of  a  four-step  reaction  mechanism.  Krewer 
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Nomenclature 

As  cross-sectional  active  electrode  area 
(1.0  X  10-3  m2) 

B  constant  depending  on  the  fuel  cell  and  its  oper¬ 

ating  state  (V) 

C  concentration  (mol  m- 3 ) 

Cfeed  feed  concentration  (mol  m-3) 

CpL  concentration  of  component  i  in  catalyst  layer 
(mol  m-3) 

Ca  double  layer  capacity  of  anode  (F  m-2) 

Cc  double  layer  capacity  of  cathode  (F  m-2) 

Cc, loss  concentration  loss  of  cathode  part  (V) 

Cf  feed  concentration  of  component  i  (mol  m-3) 

Cp, solid  heat  capacity  of  cell  (J  kg-1  K-1) 

heat  capacity  of  outlet  stream  from  anode  part 
(Jmol-1K-1) 

heat  capacity  of  outlet  stream  from  cathode  part 
(Jmol-1  K-1) 

dM  PEM  thickness  (2.0  x  10-4  m) 

Dq oh  Diffusion  coefficient  of  methanol  in  membrane 
(2.9  x  10-10  m2  s-1) 

F  Faraday  constant  (96,485  C  mol-1) 
h™  •  inlet  enthalpy  into  anode  part  ( J  mol- 1 ) 
inlet  enthalpy  into  cathode  part  (J  mol-1) 
h^f  outlet  enthalpy  from  anode  part  (J  mol-1) 
h°uf  outlet  enthalpy  from  cathode  part  (J  mol-1) 

/cell  cell  current  density  (A  m-2) 

/crossover  crossover  cell  current  density  (A  m-2) 
i\im  limiting  current  density  (A  m-2) 

kj  rate  constant  of  reaction  j  (mol  m-2  s-1) 
kLS  mass  transfer  coefficient  (ms  1 ) 
kp  hydraulic  permeability  of  membrane 
(1.57  x  10-18  m2) 

Kj  equilibrium  constant  of  reaction  j 
mso Hd  cell  mass  (kg) 

water  molecular  weight  of  water  (kg  mol-1) 
n%  mole  amount  of  anode  part  (mol) 
nc  mole  amount  of  cathode  part  (mol) 

^ch3oh  methanol  crossover  (mol  m-2  s-1) 
h™  inlet  molar  flow  rate  into  anode  part  (mol  s-1 ) 

h *n  inlet  molar  flow  rate  into  cathode  part  (mol  s-1 ) 

n&  electro-osmotic  drag  coefficient 

^ water  water  crossover  (mol  m-2  s- 1 ) 

Pa  total  pressure  of  anode  part  (Pa) 

Pc  total  pressure  of  cathode  part  (Pa) 

Po2  partial  pressure  of  oxygen  into  cathode  part  (Pa) 
Pstd  standard  pressure  of  ( 1 .0 1 3  x  105  Pa) 

Powerceii  power  generation  from  DMFC  cell  (J  s-1) 
Powercomp  compressor  power  consumption  (J  s-1) 
PowerNet  net  power  generation  from  DMFC  system 
(Js-1) 

2a,rxn  reaction  heat  of  anode  catalyst  layer  (J  s-1) 


2c,rxn  reaction  heat  of  cathode  catalyst  layer  (J  s-1) 
GMeOHXover  reaction  heat  of  methanol  oxidation  resulting 
from  methanol  crossover  (J  mol-1) 

2ioss  heat  loss  of  cell  resulting  from  convection, 
conduction,  radiation,  and  any  others  (Js-1) 
(assumed  to  be  zero  in  this  work  due  to  adiabatic 
condition) 

r  reaction  rate  (mol  m-2  s  - 1 ) 

R  universal  gas  constant  (8.314  Jmol-1  K-1) 

RME  relative  mean  error  (%) 
t  time  (s) 

Wial  start  time  of  the  operation  (s) 

/final  end  time  of  the  operation  (s) 

T  cell  temperature  (K) 

V™  inlet  flow  rate  into  the  anode  part  (m3  s-1) 

ycm  inlet  flow  rate  into  the  cathode  part  (m3  s-1) 

Vceii  cell  voltage  (V) 

yc°eu  open  circuit  voltage  of  cell  (V) 

outlet  mole  fraction  from  anode  part 
x™1  outlet  mole  fraction  from  cathode  part 

Greek  symbols 

a  charge  transfer  coefficient  (0.5) 

y  ratio  of  heat  capacity,  /.  e. ,  heat  capacity  at  constant 

pressure  divided  by  that  at  constant  volume  (1.4) 
t)sl  anode  electrode  overpotential  (V) 

y) c  cathode  electrode  overpotential  (V) 

^Fuel  feel  efficiency 

@i  surface  fraction  covered  by  component  / 
km  conductivity  of  membrane  ( 1 7  ^2  - 1  m- 1 ) 

/j  pore  fluid  viscosity  in  membrane 
(3.353  x  10-4kgm-1s-1) 

Pwater  water  density  (kg  m-3 ) 

r  mean  residence  tim  in  the  anode  or  cathode  com¬ 

partment  (s) 

Subscripts 
a  anode 

c  cathode 

feed  methanol  inlet  into  the  anode 

Pt  free  active  sites  of  platinum  catalyst 
Pt3-COH  Pt3-COH  occupied  platinum  sites 
Pt-COOH  Pt-COOH  occupied  platinum  sites 
Ru  free  active  sites  of  ruthenium  catalyst 
Ru-OH  Ru-OH  occupied  ruthenium  sites 

Superscripts 
CL  catalyst  layer 

LS  from  liquid  bulk  to  surface 

M  membrane  (PEM) 

0  at  standard  conditions 
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et  al.  [11]  presented  the  dynamic  behaviour  of  DMFCs  through 
a  dynamic  model  with  an  updated  reaction  mechanism.  Schultz 
and  Sundmacher  [12]  simulated  a  one-dimensional,  rigorous, 
mathematical  model  of  a  DMFC  with  focus  on  a  realistic  descrip¬ 
tion  of  mass  and  energy  transport,  as  well  as  the  physical 
properties  of  the  PEM  material. 

Through  the  application  of  computational  fluid  dynamics 
(CFDs),  Cheng  et  al.  [13]  derived  a  simplified  CFD  model  to  pre¬ 
dict  fuel  cell  performance  on  personal  computers  without  huge 
computation  time.  A  3-D  CFD  model  was  also  created  by  Kjeang 
et  al.  [14]  to  analyze  methanol  transport  in  a  flowing-electrolyte 
DMFC. 

Apart  from  the  above  simulation  work,  including  CFD  mod¬ 
elling  and  differential  algebraic  equation  (DAE)  modelling,  Xu 


3 

Overall  :  CH3OH(l)  +  -02(g)  ->  C02(g)  +  2H20(1) 

VCelfC  =  1  -21  V 

The  summarized  DAE  model  of  Sundmacher  et  al.  [10]  is 
listed  in  the  Table  1.  The  energy  balance  Eq.  (13),  described 
in  more  detail  in  the  Appendix  A,  includes  the  moles  of  each 
component  (liquid  and  gas  phase)  in  the  anode  and  cathode  parts, 
the  molar  flow  rates  of  the  inlet  and  outlet  streams  of  the  cell, 
the  mass  and  heat  capacity  of  the  solid  state  (the  single-cell 
stack),  the  reaction  heat  of  methanol  oxidation  due  to  methanol 
crossover,  and  the  electrochemical  reaction  heats  of  anode  and 
cathode  catalyst  layers. 


d  T 
d  t 


[e, 


(hfj  -  h°f)  +  E,n"’-van,(/'an,  -  KUb  -  2c,rxn  -  2a,rxn  -  QMeOHXover«CH3OHAS  “  Powercell  -  Q\, 


WIsolidCpsolid  +  racEi(<?Cpc!;)  +  «aE,(*a“tCpa!7 


(13) 


et  al.  [15]  performed  optimization  of  the  isothermal  DMFC  DAE 
model  based  on  the  study  by  Sundmacher  et  al.  [10]. 

In  addition  to  the  active  research  on  DMFC  modelling  and 
the  strong  influence  of  membrane-electrode  assembly  (MEA) 
temperature  on  DMFC  performance,  an  optimization  study  of 
a  non-isothermal  DMFC,  which  involves  the  flow  rates  and 
the  reaction  enthalpies  of  inlet  and  outlet  streams,  is  neces¬ 
sary  to  determine  an  optimum-operating  condition.  Thus,  the 
work  reported  here  presents  the  development  of  a  non-isothermal 
dynamic  DMFC  model  by  adding  an  energy  balance  to  the 
isothermal  model  of  Sundmacher  et  al.  [10],  together  with 
dynamic  simulations  and  optimization  of  the  model.  The  next 
section  introduces  a  non-isothermal  differential  algebraic  equa¬ 
tion  model,  updated  from  the  isothermal  model  of  Sundmacher 
et  al.  [10],  and  parameter  estimation  results.  The  dynamic 
behaviour  of  DMFC  operation,  i.e.,  the  effects  of  operating  vari¬ 
ables,  and  a  dynamic  optimization  are  discussed  in  the  following 
sections. 


Based  on  the  energy  balance  the  outlet  temperature  from  the 
DMFC  can  be  predicted,  under  given  conditions  of  the  inlet 
temperatures,  molar  flow  rates  and  feed  concentration,  by  adopt¬ 
ing  some  of  the  assumptions  made  by  Sundmacher  et  al.  [10], 
namely: 

•  the  compartment  channels  and  catalyst  layers  are  treated  as 
continuous  stirred  tank  reactors  (CSTR),  i.e.,  a  lumped  MEA 
model 

•  ohmic  drops  in  current-collectors  and  electric  connections  are 
negligible 

•  mass-transport  resistances  in  the  catalyst  layers  are  negli¬ 
gible  because  the  catalyst  layers  are  very  thin  (i.e.,  10p,m) 
compared  with  that  of  the  PEM  (i.e.,  200  p,m) 

•  mass-transport  coefficients  of  the  components  (water  and 
methanol)  in  the  diffusion  layer  are  equal 

•  the  electrochemical  oxidation  of  methanol  at  the  anode  fol¬ 
lows  the  reaction  mechanism  given  in  Table  1  with  the  first 
step,  Eq.  (1),  as  a  rate-determining  step. 


2.  Non-isothermal  dynamic  model  of  DMFC 

2.1.  Differential  algebraic  equation  (DAE)  model  of  DMFC 

The  basic  structure  of  a  DMFC  reaction  scheme  is  shown 
in  Fig.  1,  and  the  following  total  electrochemical  reactions  take 
place  on  the  catalyst  layers  of  the  MEA. 

Anode  :  CH3OH(l)  +  H20(1)^*C02  +  6H+  +  6e“ 

y#,25  °C  —002V 
VCH30ll±>C02  —  U'UZ  v 

Cathode  :  -02(g)  +  6H+  +  6e“-2>3H20(l) 

25°c  —  1  ?3V 

V02**H20  —  A’ZJ  V 


Also,  the  additional  assumptions  are  made: 

•  there  is  no  heat  loss  in  the  stack  (adiabatic  condition),  that  is, 

2loss  —  0 

•  the  electro-osmotic  drag  coefficient  (nf)  is  a  function  of 
the  temperature  (^a  = —5.77 +  0.027  7);  this  linear  equa¬ 
tion  is  obtained  by  fitting  the  data  of  Ren  and  Gottesfeld 
[16] 

•  the  inlet  flow  rate  into  the  anode  or  the  cathode  is  the  same 
as  the  respective  outlet  flow  rate  based  on  the  lumped  MEA 
model  assumption. 

Water  crossover  is  derived  from  the  paper  of  Lu  et  al.  [17], 
i.e., 

=  <»d  +  0.5)  (‘-f)  -  V.  -  (14) 

\  r  y  [1  ^Wwater 
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Membrane 

Electrode 

Assembly 


Fig.  1.  Structure  and  reaction  scheme  of  DMFC. 


Table  1 

DAE  model  taken  from  Sundmacher  et  al.  [10] 


Electrode  reaction  rate  expressions 
Anode 


r\  —  k\  exp  ^  ^  l^ptCcHsOH  fr  exP  ^  ^  6>Pt3-C0H| 

r2  =  k2  exp  j^Ru  -  Y2  exp  (_ /^a)  0Ru-°H} 

r3  =  k2  1 0Pt3-COH  ^Ru-OH  _  —  6)pt-COOH@Pt0Ru| 

|  0pt-COOH@Ru-OH  -  —  Cqq2  0pt0Ru  | 


r4  =  k4  0pt— COOH0RU-OH  “ 

Cathode 

fa5F 

r5  =k5  exp  ( 

Mass  balances  of  anode  compartment 
dCcH3OH  _  ^ch3oh  _  cch3oh  _  khSAs 
d  t  Ta  Va 

dCco2  _  ^co2  —  ^c°2  khSAs 

d  t  ra  va 

Mass  balances  of  catalyst  layer 

dCCH3OH  _  kLSAS  rO.  x  "  M  - 

—  ,/CL  'lch3oh  <-ch3oh-*  1/CL  "ch3oh  ,/CL ri 

va  va  va 


(Cch3OH  CcH3Oh) 


(Cco2  C^q2  ) 
As 


As 


dC 


d  t 

-CL 


C02  _  kLSAs  CL  As 

(Cc02  —  Cco  )  + 


'2  -C(V-r  CLri 
V  a 


dt  V' 

Charge  balances 
dr] a  _  icell  ~  dFri 
dt  ~  Ca 

dr]c  _  -icell  ~  0F(r$  +  n^CH30H) 

~di  ~  C~c 

The  methanol  crossover 


/7m 

"ch3oh  — 


DCH3OHC£cH3OH  \  /  exp(Pg)  \ 

dM  )  6  \exp(Pe)  —  1  / 


here,  Pe  =  v- 


D: 


(M 

ch3oh 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(10) 

(ID 

(12) 
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where  kv  is  the  hydraulic  permeability  of  a  membrane;  p  is  the 
pore  fluid  viscosity  in  the  membrane;  pWater  is  the  water  density; 
Mwwater  is  the  molecular  weight  of  water. 

Since  the  cathode  overpotential  does  not  contain  concentra¬ 
tion  loss  whereas  the  anode  overpotential  does,  the  concentration 
loss  at  the  cathode  is  given  by  [18]: 

Cc>loss  =  51n(l-^!l)  (15) 

\  flim  / 

where  B  is  a  constant  depending  on  the  fuel  cell  and  the  operating 
state  [18]. 

The  overall  cell  voltage  is  modified  from  that  of  the  paper 

[10], 

Ken  =  vc°ell  -  rla  +  >1c  -  ^cell  +  B  In  (l  -  — )  (16) 

K  \  Aim  / 

2.2.  Optimization  strategy 

Based  on  the  DAE  model  in  Section  2.1,  parameter  estima¬ 
tion,  dynamic  simulation  and  optimization  were  performed  by 
using  the  gPROMS  modelling  tool  [  19] .  For  the  dynamic  simula¬ 
tion,  a  DASOLV  code,  one  of  the  standard  differential  algebraic 
solvers  in  the  gPROMS,  was  employed  with  a  dynamic  opti¬ 
mization  solver,  DYNOPT,  for  the  optimization  and  parameter 
estimation. 

The  DASOLV  solver,  based  on  a  variable  time  step  and 
variable  order  backward  differentiation  formulae  (BDF),  is 
applicable  for  a  wide  range  of  problems  due  to  its  efficient  treat¬ 
ment  of  large  and  sparse  systems  of  equations  that  include  the 
variable  within  the  lower  and  upper  bounds  except  highly  oscilla¬ 
tory  behaviour  models  with  frequent  discontinuities,  It  provides 
strong  advantages  in  automatically  adjusting  each  time  step  and 
in  yielding  efficient  finite  difference  approximations. 

A  CVP_SS  solver  of  the  DYNOPT,  adopted  for  the  opti¬ 
mization,  is  based  on  a  control  vector  parameterization  (CVP) 
approach  assuming  that  the  time-varying  control  variables  are 
piecewise  constant  (or  piecewise  linear)  functions  of  time  over 
a  specified  number  of  control  intervals. 

The  CVP_SS  solver  employs  a  single- shooting  dynamic  opti¬ 
mization  algorithm,  that: 


Table  2 


Reference  values  of  kinetic  parameters 


Parameters 

Reference  values 

11 

-j* 

1.0  x  10“7 

kl  —  ^2,  ref 

1 

&3  =  ^3, ref 

0.01 

^4  =  ^4,  ref 

1 

1 

II 

>/-> 

7.0  x  10-5 

Kl  —  K\  ref 

1 

k2  =  k2M 

0.02 

K3=  K3,  ref 

1 

K4  =  £4,  ref 

0.02 

2.3.  Parameter  estimation 

The  parameters  to  be  estimated  are  selected  from  the  chemi¬ 
cal  reaction  parameters,  i. e. ,  k\ ,  &2,  &3 ,  &4,  &5 ,  K\ ,  ,  and  K4 

in  the  Eqs.  (1)— (5)  of  Table  1.  The  three  parameters  (&2,  &3,  and 
k/f)  are  screened  in  advance  because  they  have  no  influence  on 
the  DMFC  behaviour  in  the  current  model,  i.e.,  r2  =  r3  =  r\  ~  0 
based  on  the  chemical  equilibrium  assumption  and  allow  cal¬ 
culation  of  the  surface  fraction  of  catalyst  areas  [10].  Then,  the 
sensitivities  of  the  remaining  parameters  (k\,  £5,  K\,  K2 ,  ^3, 
and  K4)  to  the  cell  voltage  are  analyzed  to  choose  those  for  esti¬ 
mation  by  using  the  reference  values  in  Table  2.  In  Table  3,  Vj 
is  the  voltage  when  the  chemical  reaction  parameter  is  /,  with 
J  €=  {^4,ref>  &5, ref?  ^T,ref>  •  •  -v^4, ref}-  Vjx  0.5  Vj  x  1.5  the 

voltages  when  the  corresponding  parameter  value  is  multiplied 
by  0.5  and  1.5,  respectively. 

Since  the  cell  voltage  is  mainly  affected  by  four  kinetic 
parameters  (k\ ,  £5 ,  K2,  and  K4),  as  shown  in  Table  3,  these  param¬ 
eters  are  finally  selected  as  decision  variables  for  the  estimation. 

The  constant  ( B )  and  limiting  current  density  (z)im)  in  the 
cathode  concentration  loss  term  are  also  treated  as  parameters  for 
estimation.  The  optimization  model  for  the  parameter  estimation 
is  formulated  as: 


(i)  chooses  the  duration  of  each  control  interval  and  the  values 
of  the  control  variables 

(ii)  solves  a  dynamic  system  model  over  the  entire  time  horizon 
by  starting  from  the  initial  time  point 

(iii)  determines  the  values  of  the  optimized  objective  function 
to  satisfy  the  constraints 

(iv)  repeats  the  procedure  until  convergence  to  the  optimum  is 
reached. 


The  non-linear  programming  (NLP)  solver  of  the  CVP_SS 
is  a  SRQPD  solver  which  employs  a  sequential  quadratic  pro¬ 
gramming  (SQP)  approach. 


Table  3 

Sensitivity  analysis  of  kinetic  parameters  to  voltage  (zceii  =  2000  A  m-2) 


Parameters 


Voltage  deviation  [%]  when  /=ki,ref, 
•  •  •?  ^5,  ref?  A 1  ^ef ?  .  .  Or  ^4;ref 


V/xO.5  — V/ 

Vj 


x  100 


Vjxijs-Vj 

Vj 


X  100 


J=k  1 

2.58 

1.39 

J=k5 

10.73 

6.28 

J=KX 

-0.00 

-0.00 

j=k2 

4.18 

2.41 

j=k3 

-0.00 

-0.00 

J=Ka 

4.18 

2.41 

76 


D.  Ko  et  al.  /  Journal  of  Power  Sources  180  (2008)  71-83 


Table  4 

Initial  guess  values  for  parameter  estimation 


Parameters 


Reference  values 


h 

h 

K2 

k4 

B 

*lim 


1.0  x  10  7  (molm  2  s  ') 
7.0  x  10-5  (molm-2  s-1) 


0.02 
0.02 
0.01  (V) 
7000  (Am-2) 


Subject  to: 


1(T4V  <B<  IV 
6000  Am-2  <  i\{m  <  10,  000  Am-2 
1 .0  x  10-9  mol  m-2  s_1  <  k\  <  1 .0  x  10-5  mol  m-2  s_1 
1.0  x  10-10molm-2s-1  <  ks  <  1.0  x  10-4  molm-2  s_1 
1.0  x  10“4  <  K2  <  5.0 
1.0  x  10“4  <  K4  <  5.0 


Model  Eqs.  (1)— (16),  optimization  tolerance  =  10-6. 

The  initial  guess  values  for  this  parameter  estimation  are 
shown  in  Table  4. 

In  Eq.  (17) /is  an  objective  function  to  be  minimized;  B  is  a 
constant  that  depends  on  the  fuel  cell  and  its  operating  state  [18]; 
him  is  the  limiting  current  density;  /^ell  and  ilce[l  are  the  current 
densities  of  the  k-th  and  Z-th  cell,  respectively;  and 

VCEJ O'Neil)  are  calculated  and  measured  cell  voltage,  respec¬ 
tively,  when  the  curent  density  is  /£el];  7?aal(4eii)  an(*  ^aXP(4ell) 
are  the  calculated  and  measured  overpotential  of  the  anode, 
respectively,  when  the  current  density  is  i4  The  measured 
number  of  data  points  for  cell  voltage  is  thirty  and  that  for  anode 
overpotential  is  three. 

The  estimated  values  of  the  parameters  are  listed  in  Table  5. 
The  relative  mean  errors  (RMEs)  of  the  cell  voltage  and  anode 


Table  5 

Estimated  values  for  parameters 


Parameters  Reference  values 


k\ 

1.76  x  1  cr7  [molm- 

h 

1.00  x  10“4  [mol  m" 

k2 

3.45  x  10"2 

k4 

3.45  x  10“2 

B 

1.00  x  10"4  [V] 

*lim 

7034.84  [Am'2] 

Total  CPU  time  for  calculation:  552.94  s  on  Pentium  4  with  3  GHz  CPU  and 
512  MB  RAM,  number  of  NLP  iterations:  8. 


Fig.  2.  Parameter  estimation  results:  cell  voltage  and  anode  overpotential  vs. 
current  density  (cell  temperature  =  343.15  K). 


overpotential  are  calculated  by  Eqs.  (18)  and  (19). 


RMEVcell  = 


1 

30 


vCal //*  \  _  VExp(ik 
’cell  *7  cell!  vcell  UCelU 


yExP  (jk  \ 

Kcell  UcelU 


3.51% 

(18) 


RME^a 


^cefl^cell)  ^cell  ^cell) 


Exp,./ 


4.40%  (19) 


Because  the  RMEs  are  reasonably  small  and  the  simulated 
data  are  in  good  agreement  with  the  measured  data,  as  shown 
in  Fig.  2,  the  estimated  values  of  the  model  parameters  describe 
the  steady- state  performance  of  the  DMFC  well. 


3.  Dynamic  behaviour  of  DMFC 

The  effect  of  the  operating  variables  (feed  concentration, 
anode  flow  rate,  and  cathode  flow  rate)  on  the  temperature, 
crossovers  and  voltage  was  checked  using  the  obtained  param¬ 
eter  values,  when  the  current  density  is  2000  Am-2.  The  feed 
concentration  has  an  influence  on  the  cell  voltage,  temperature 
and  crossovers  of  methanol  and  water,  as  shown  in  Figs.  3-6. 
The  temperature  increases  as  the  feed  concentration  increases 
(Fig.  3)  because  the  higher  feed  concentration  leads  to  higher 
methanol  crossover  (Fig.  4).  The  methanol  crossover  also  influ¬ 
ences  the  methanol  oxidation  reaction,  which  is  exothermic.  The 
feed  concentration  affects  the  water  crossover  too,  as  shown  in 
Fig.  5,  due  to  its  effect  on  the  system  temperature.  The  cell  volt¬ 
age  decreases  as  the  feed  concentration  increases  at  2000  A  m-2 
(Fig.  6)  due  to  the  strong  influence  of  crossover  on  the  cathodic 
overpotential. 

As  shown  in  Fig.  7,  higher  flow  rate  of  feed  lowers  the  cell 
temperature  at  constant  feed  concentration  (1  M).  This  means 
the  feed  acts  as  a  coolant.  Methanol  crossover  decreases  as  the 
flow  rate  decreases  since  the  amount  of  methanol  supplied  is 
reduced  (Fig.  8).  Due  to  the  effect  of  temperature  on  the  water 
drag  coefficient,  water  crossover  increases  as  the  flow  rate  of 
feed  decreases  (Fig.  9).  The  cell  voltage  changes  proportionally 
to  the  flow  rate  (Fig.  10). 
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Fig.  3.  Temperature  profile  vs.  feed  concentration  at  non-isothermal  conditions  Fig.  5.  Water  crossover  profile  vs.  feed  concentration  (initial  tempera- 


(initial  temperature  =  25  °C,  current  density  =  2000  A  m 
ature  =  333.15  K,  cathode  inlet  temperature  =  298.15  K). 


■,  anode  inlet  temper- 


ture  =  25  °C,  current  density  =  2000  A  m 
cathode  inlet  temperature  =  298.15  K). 


■,  anode  inlet  temperature  =  333. 15  K, 


As  the  air  flow  rate  increases,  the  temperature  decreases 
because  air  also  affects  as  a  coolant  (Fig.  11).  The  crossovers  of 
methanol  and  water  decrease  when  the  air  flow  rate  increases, 
as  demonstrated  in  Figs.  12  and  13.  The  data  in  Fig.  14  show 
that  a  higher  air  flow  rate  results  in  a  higher  cell  voltage. 


conditions  such  as  inlet  flow  rates  into  the  anode  and  the  cath¬ 
ode,  as  well  as  the  feed  concentration.  Fig.  15  shows  the  assumed 
change  in  power  density  load.  The  optimization  model  is  formu¬ 
lated  as  follows: 


4.  Dynamic  optimization  of  DMFC 

A  dynamic  optimization  of  the  DMFC  was  performed  to 
achieve  the  desired  fuel  efficiency,  80%,  and  to  satisfy  the 
varying  power  load.  The  decision  variables  are  the  operating 
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Fig.  4.  Methanol  crossover  profile  vs.  feed  concentration  (initial  tempera¬ 
ture  =  25  °C,  current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333.15  K, 
cathode  inlet  temperature  =  298. 15  K). 
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Fig.  6.  Cell  voltage  profile  vs.  feed  concentration  (initial  temperature  =  25  °C, 
current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333.15  K,  cathode  inlet 
temperature  =  298.15  K). 
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Fig.  7.  Temperature  profile  vs.  feed  flow  rate  at  non-isothermal  condition 
(initial  temperature  =  25  °C,  current  density  =  2000  Am-2,  anode  inlet  temper¬ 
ature  =  333.15  K,  cathode  inlet  temperature  =  298.15  K). 
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Fig.  9.  Water  crossover  profile  vs.  feed  flow  rate  (initial  temperature  =  25  °C, 
current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333.15  K,  cathode  inlet 
temperature  =  298. 15  K). 


<  Vln  <  1.0  X  10 


1  5  m3  s  1 


Subject  to: 

1.0  x  10-8 
1.0  x  10-7  <  Vcm  <  1.0  x  10_4m3s_1 
900  <  Cfeed  <  1100  mol  m-3 

Model  Eqs.  (1)— (26),  optimization  tolerance  =  10_3.Here, 

/  crossover  =  (6F)(«ch3Oh) 

/cell 


where  rj^uh  c^lcFilated  fuel  efficiency;  rj^uei  ave 

average  calculated  fuel  efficiency;  ^pug}red  is  the  desired  fuel 
efficiency;  /crossover  is  the  crossover  current  density;  ffinai  is  the 
end  time  of  the  operation;  initial  is  the  start  time  of  the  operation 
(=0  s).  The  power  equations  are  described  by  Eqs.  (24)— (26). 
The  compressor  power  consumption  (Powercomp),  assumed  as 
an  adiabatic  condition,  is  given  by  Eq.  (24)  [20]. 


„Cal 


is  the 


Cal  __ 
^Fuel 


^Fuefave 


/  cell  +  ^crossover 

fiinal  initial 


Methanol  Crossover  [mol/mA2/sec]  (left) 


(21) 

(22) 

(23) 


- Feed  Flowrate  [mA3/sec]  (right) 

Fig.  8.  Methanol  crossover  profile  vs.  feed  flow  rate  (initial  temperature  =  25  °C, 
current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333. 15  K,  cathode  inlet 
temperature  =  298. 15  K). 


Powercomp  — 


371  rcm]/ycm 

(V  —  l)^comp,eff 


-) 
v  ^STD  ) 


(1— (i/y)) 


—  1 


(24) 


where  y  is  the  ratio  of  heat  capacity,  i.e.,  the  heat  capacity  at 
constant  pressure  divided  by  that  at  constant  volume  (=1.4); 
Tcm  is  the  inlet  temperature  into  the  cathode  ;  V™  is  the  inlet 
flow  rate  into  the  cathode;  ^Comp,eff  is  compressor  efficiency 
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Fig.  10.  Cell  voltage  profile  vs.  feed  flow  rate  (initial  temperature  =  25  °C,  cur¬ 
rent  density  =  2000 Am-2,  anode  inlet  temperature  =  333.15 K,  cathode  inlet 
temperature  =  298. 15  K). 
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Temperature  [K]  (left)  - Air  Flowrate  [mA3/sec]  (right) 


Fig.  11.  Temperature  profile  vs.  air  flow  rate  at  non-isothermal  condition 
(initial  temperature  =  25  °C,  current  density  =  2000  Am-2,  anode  inlet  temper¬ 
ature  =  333.15  K,  cathode  inlet  temperature  =  298.15  K). 


(=80%);  Pc  is  total  pressure  at  the  cathode;  Pstd  is  the  standard 
pressure. 

The  power  generation  of  the  DMFC  cell  (Powerceii)  is  calcu¬ 
lated  by: 

Powercen  =  /Ceii  x  Vceii  x  (25) 

Thus,  the  net  power  generation  from  the  DMFC  system 
(PowerNet)  is  obtained  by  Eq.  (26),  and  it  must  be  the  same 
as  the  required  power  load. 

PowerNet  =  Powerceii  -  PowerCOmp  (26) 

The  initial  guess  values  for  the  optimization  are  listed  in 
Table  6.  The  optimum  dynamic  operating  conditions  of  the 
feed  and  air  flow  rates,  and  feed  concentration  are  given  in 
Figs.  16-18,  respectively.  The  data  show  that  the  decision  vari- 
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Fig.  12.  Methanol  crossover  profile  vs.  air  flow  rate  (initial  temperature  =  25  °C, 
current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333. 15  K,  cathode  inlet 
temperature  =  298. 15  K). 
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Fig.  13.  Water  crossover  profile  vs.  air  flow  rate  (initial  temperature  =  25  °C, 
current  density  =  2000  A  m-2,  anode  inlet  temperature  =  333. 15  K,  cathode  inlet 
temperature  =  298.15  K). 
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Fig.  14.  Cell  voltage  profile  vs.  air  flow  rate  (initial  temperature  =  25  °C,  cur¬ 
rent  density  =  2000 Am-2,  anode  inlet  temperature  =  333.15 K,  cathode  inlet 
temperature  =  29 8 . 1 5  K) . 
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Fig.  15.  Assumed  power  load  in  non-isothermal  DMFC  single-cell  model. 


80 


D.  Ko  et  al.  /  Journal  of  Power  Sources  180  (2008)  71-83 


Table  6 

Initial  guesses  for  decision  variables  in  optimization 


1.00E+03  - 


Variables 

Initial  guess  value 

9.99E+02 

Feed  flow  rate  (m3  s-1) 

1.7  x  1CT7 

Air  flow  rate  (m3  s-1 ) 

2.0  x  10“5 

JT  9.94E+02 
e 

Feed  concentration  (molm-3) 

1000 

c 

o 

i 

O' 
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O 


200 
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Time  [sec] 
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Fig.  16.  Feed  flow  rate  in  non-isothermal  DMFC  model  before  and  after  opti¬ 
mization. 

ables  (flow  rate  of  feed  and  air,  and  feed  concentration)  change 
proportionally  to  the  power  density  load.  That  is  the  optimum 
values  of  the  operating  variables  increase  when  the  power  load 
density  increases.  The  optimum  feed  concentration  is  lower  than 
the  initial  guess  value  because  the  methanol  crossover  is  reduced 
with  the  low  feed  concentration  and  leads  to  high  fuel  efficiency. 
The  optimum  air  flow  rate  increases  compared  with  the  value 
before  the  optimization,  because  a  higher  rate  results  in  a  higher 
voltage. 
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Fig.  18.  Feed  concentration  in  non-isothermal  DMFC  model  before  and  after 
optimization. 
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Fig.  17.  Air  flow  rate  in  non-isothermal  DMFC  model  before  and  after  opti¬ 
mization. 


Fig.  19.  Fuel  efficiency  before  and  after  optimization  (desired  fuel  effi¬ 
ciency  =0.80). 


Judging  from  the  results  of  the  dynamic  optimization  and 
simulations,  we  find  that  methanol  crossover  is  more  affected 
by  feed  concentration  than  by  the  flow  rates  of  the  feed  and  air. 

In  this  study,  the  optimum  objective  function  value  is 
3.23  x  10-5  and  the  average  fuel  efficiency  is  79.99%,  which  is 
extremely  close  to  the  desired  fuel  efficiency  (=80%)  as  shown 
in  the  Fig.  19.  Given  these  optimum  conditions,  the  average  fuel 
efficiency  is  enhanced  from  71.78  to  79.99%.  The  total  CPU 
time  for  the  optimization  is  1213.98  s  on  Pentium  4  with  3  GHz 
CPU  and  512  MB  RAM  and  the  number  of  NLP  iterations  is  8. 

5.  Conclusions 

A  non-isothermal  dynamic  DMFC  model  is  developed  and 
incorporates  electrochemical  reaction  kinetics  and  crossover 
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phenomena.  This  study  is  a  foundation  work  to  predict  accu¬ 
rate  stack  temperature  and  dynamic  behaviour  for  the  efficient 
design  of  a  DMFC  stack  and  to  find  the  optimum-operating 
conditions.  For  the  further  analysis  and  optimization  of  the 
non-isothermal  dynamic  fuel  cell  behaviour,  reliable  heat  loss 
(convection,  conduction,  and  radiation)  needs  to  be  considered 
within  a  non-adiabatic  and  non-isothermal  model.  Thus,  an 
accurate  temperature  of  a  DMFC  stack  can  be  determined  by 
including  the  exact  heat  loss  calculation  through  an  extended 
model.  An  accurate  temperature  prediction  should  be  very  help¬ 
ful  in  reducing  the  cost  of  the  optimum  design  and  operation. 
The  current  basic  model  is  adopted  from  the  isothermal  DMFC 
model  of  Sundmacher  et  al.  [10]  by  adding  the  energy  balance 
equation  to  predict  the  influence  of  temperature  on  the  concentra¬ 
tion  loss  of  the  cathode  catalyst  layer  and  on  water  crossover.  To 
fit  the  non-isothermal  DMFC  MEA  model  to  a  practical  MEA, 
model  calibration  through  the  optimization  technique  is  carried 
out  by  estimating  the  chemical  reaction  constants  (k\  ,  £5 ,  K2,  and 
K4)  and  the  parameters  (B  and  fiim)  of  the  cathode  concentration 
loss  term.  The  I-V  and  I-r /a  points  of  the  calibrated  model  show 
good  agreement  with  the  experimental  data.  The  dynamic  sim¬ 
ulation  and  optimization  of  the  calibrated  model  are  carried  out 
with  a  gPROMS  modelling  tool  that  includes  the  SRQPD  algo¬ 
rithm  for  an  optimization  and  a  DASOLV  code  for  a  dynamic 
simulation.  The  present  model  describes  the  dynamic  DMFC 
behaviour  (voltage,  crossover  of  methanol  and  water,  temper¬ 
ature)  successfully  according  to  step  changes  of  the  operating 
variables  (feed  concentration,  inlet  flow  rate  into  anode  part, 
inlet  flow  rate  into  cathode  part).  The  simulation  results  show 
that  the  feed  concentration  is  a  very  important  operating  vari¬ 
able  -  because  of  its  strong  influence  on  the  dynamic  behaviour 
of  the  DMFC  (temperature,  crossover  -  voltage)  compared  with 
the  flow  rates  of  the  feed  and  air.  Also,  the  optimum-operating 
strategy  can  easily  be  established  to  achieve  the  desired  fuel  effi¬ 
ciency  and  to  satisfy  the  varying  power  load  through  the  dynamic 
optimization  work.  In  conclusion,  the  simulation  and  optimiza¬ 
tion  can  efficiently  provide  an  insight  into  the  DMFC  during 
dynamic  operation  and  suggest  optimum- operating  strategies 
without  the  need  for  experimental  data. 

Acknowledgements 

The  authors  are  grateful  to  Dr.  Heetak  Kim  and  Dr.  Junyoung 
Park  at  the  fuel  cell  group  of  Samsung  SDI  for  providing  exper¬ 
imental  data  and  helpful  discussion  on  the  modelling  results. 


Appendix  A.  Description  of  energy  balance 

The  energy  balance  Eq.  (13)  consists  of  heat  inlet  and  outlet, 
reaction  heats,  heat  loss,  and  methanol  oxidation  heat  resulting 
from  the  crossover.  In  this  section,  the  equations  included  in  the 
energy  balance  are  described.  The  heat  capacities  of  the  cathode 
and  the  anode  ( Cpc  and  Cp a,  respectively)  are  given  by  [21]: 

=  Ai  +  B,T  +  C,T2  +  D,T~2  (Al) 

R 


=  A,  +  BjT  +  Ci  T2  +  Di  T~2  (A2) 

R 

where,  the  A/,  B[ ,  Q,  and  Dl  are  the  coefficients  of  heat  capacity 
as  listed  in  Table  Al. 

The  reaction  heats  of  the  cathode  and  the  anode  (<2c,rxn  and 
Ga,rxn,  respectively)  are  calculated  by: 

Gc,rxn  —  ^  ,i  r5  A  )  (A3) 


<2a,rxn  =  iAS)  (A4) 


The  enthalpies  of  the  inlet  and  the  outlet  streams  are  given 

by  [21], 


h\ "•  =  [  C  R  (  C]^i- )  dT 

JrKf= 298. 15K  L  l  R  )  . 
cTl 


/ 

J  rref=298.15  K 


+  BiT  +  CjT2  +  DjT~2}]dT 


_  r/Cathode  ,  r> 
—  ^298  ,i  +  K 


Aj  +  Bi 


+  rref 


+jl(Tc'  +  rrey  -  (rcinrref)}  +  — 


Dt 


(T™Tref) 


x(72n  -  Ti-ef) 


(A5) 


Table  Al 

Coefficients  of  heat  capacity  for  enthalpy  calculations  [20] 


Component  i 

A 

fix  to3 

Cx  106 

D  x  10“5 

Temperature  range  (K) 

Anode 

CH3OH 

1=1 

13.431 

-51.28 

131.13 

0 

273-373 

h2o 

i  =  2 

8.712 

1.25 

-0.18 

0 

273-373 

C02 

i  =  3 

5.457 

1.045 

0 

-1.157 

298-2000 

Cathode 

O2 

i  =  l 

3.639 

0.506 

0 

-0.227 

298-2000 

n2 

i  =  2 

3.280 

0.593 

0 

0.040 

298-2000 

h2o 

i  =  3 

8.712 

1.25 

-0.18 

0 

273-373 
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Ki  = 


/rref=298.15K 

rTm 

/rref=298.15K 


R 


d  T 


[R{Ai  +  BJ  +  CiT 2  +  DiT~2}]dT 


Tj Anode  >  z? 
M298,i  +  K 


Ai  +  Bj 


T™  +  Af 


+  y  (CC  +  TrJ  -  (Tj"  Af)}  +  7-^- 
3  Uanrref) 


Table  A2 

Joules  per  mole  of  substance  formed  at  298  K  (standard  heat  of  generation)  [20] 
Component  i  #298,;  State 


Anode 


CH3OH 

i  =  1 

h2o 

i  =  2 

co2 

i  =  3 

Cathode 

O2 

i  =  l 

n2 

i  =  2 

h2o 

i  =  3 

-238,660 

Liquid 

-285,830 

Liquid 

-393,509 

Gas 

0 

Gas 

0 

Gas 

-393,509 

Liquid 

x(r“  -  Tref) 


(A6) 


yOUt 


Tref= 298. 15  K 

yOUt 


Kl  =  [ 

r^r 

Jtm= 298. 15  K 

_  zjCathode  i  p 

—  ^298,  i  +  ^ 

c, 


7? 


dr 


A,-  +  Bi 


7?ut  +  rref 


A; 


+y((T0Ut  +  Ta ft  -  (TcoatTre f)}  +  Tk[) 

x(7;out  -  Af) 


(A7) 


Cp.mh.MeOH  Oxid  =  3 ,9915+(55 .584  x  1(T3)  ^r  +  298-15 

R 

(-131.49  x  10“9)  , 

+  - - -{(T  +  298.15)2 


2k/T| 

(-0.8165  x  105) 

298.15  7}  + 

298.15  T 

(A  10) 

TjCell,  total  tj  Anode,  total  ,  /^Cathode,  total 

**298  ~  ^298  "h  ^298 

(All) 

-| 

tj  Anode,  total  tj  Anode  tj  Anode  zj  Anode 

**298  —  ^298,/=3  “  n29%,i=2  n29%,i=\ 

(A12) 

TjCathode, total  o  zjCathode  zjCathode  i  r  tj  _  i  Cathode 

H298  —  d#298,i=3  M29%,i=2  ~  ES#298,i  1 


_  o  r/^aim 

—  2>/7oqq 


Cathode 


298,  i=3 


(A  13) 


K1=  I7**  [/?{^i}ldr 

J  7’ref=298. 1 5  K  L  l  *  J_ 

^yOUt 


L 


[7?{A,  +  s(r  +  c,r2  +  Ar_2}]dr 


Tref =298. 15  K 

M  + 


zj Anode  i  p 
**298, i  +  K 


raout  +  A  f 


+y{(raout  +  Af)2  -  (raoutAf)}  + 
x(raout  -  Af) 


(A8) 


Here,  rref  =  298.15  K  and  the  values  of  the  parameters  (A/, 
Bi,  Ci ,  and  £>;)  are  shown  in  Table  Al. 

The  outlet  temperatures  from  the  cathode  and  the  anode  (Tcout 
and  T°ut,  respectively)  are  assumed  to  be  the  same  as  the  cell 
temperature  (T).  Tcin  and  T[n  are  the  inlet  temperatures  into  the 
cathode  and  the  anode. 

The  methanol  oxidation  heat  (GMeOHXover)  is  calculated  by 
the  Eq.  (A9),  and  the  methanol  heat  capacity  (CAmh,MeOH  Oxid) 
and  total  enthalpy  change  (/f^e811,total)  by  the  Eqs.  (A10)  and 
(All),  respectively  [21]. 


GMeOHXover  =  +  Cp,mh,MeOH  Oxid(^  ~  ?ref)  (A9) 


Each  value  of  #298,/  is  shown  in  Table  A2. 
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